library(repro)
# load packages from yaml header
automate_load_packages()
# include external scripts
automate_load_scripts()
# load data
data_ambroise_full <- automate_load_data(data_ambroise_full, read_sav)
sessio = sessionInfo(); #opts_chunk$set(echo = F, message=F, warning=F) # set echo F for all
This file was automatically created via the Repro package (version 0.1.0) using R version 4.0.1 (2020-06-06)
options(scipen = 666, warn=-1, contrasts=c("contr.sum","contr.poly"), knitr.kable.NA = '') #remove scientific notation # remove warnings #set contrasts to sum !
set.seed(666) #set random seed
# panderOptions('knitr.auto.asis', FALSE) #remove auto styling
# Look at R/clean.R (listed in the YAML) which does all the preprocessing for more info
# If you are unsure weather or not you have `git` `make` & `docker`.
# check_git()
# check_make()
# check_docker()
Blabla
base[c("Age", "Gender", "Profession")] %>% tbl_summary(statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} / {N} ({p}%)"),) %>% modify_caption("**Table 1. Demographics **") %>%
bold_labels()
| Characteristic | N = 3031 |
|---|---|
| Age | 21.7 (4.4) |
| Gender | |
| Homme | 56 / 303 (18%) |
| Femme | 244 / 303 (81%) |
| Autre/NA | 3 / 303 (1.0%) |
| Profession | |
| Etudiant.e | 298 / 303 (98%) |
| Actif.ve | 3 / 303 (1.0%) |
| Les deux | 2 / 303 (0.7%) |
|
1
Mean (SD); n / N (%)
|
|
base_clean[c("Decision Mode")] %>% tbl_summary(statistic = list(all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} / {N} ({p}%)"),) %>% modify_caption("**Table 2. Decision Mode **") %>%
bold_labels()
| Characteristic | N = 2961 |
|---|---|
| Decision Mode | |
| Affectif | 32 / 296 (11%) |
| Les deux | 202 / 296 (68%) |
| Cognitif | 62 / 296 (21%) |
|
1
n / N (%)
|
|
\(~\) \(~\)
\(~\) \(~\)
base_price = base_clean #no huge outliers, see appendix
final_price = lm(Price ~ Priming*Product + Age + Decision_mode + Priming*Product + Product*negatif_affect + Product*positive_affect, data = base_price)
apa_lm <- apa_print(anova(final_price))
apa_table(apa_lm$table, caption = "Anova table for Price.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .001 | [.000, .017] | 0.36 | 1 | 285 | .552 |
| Product | .009 | [.000, .037] | 2.72 | 1 | 285 | .100 |
| Age | .007 | [.000, .032] | 2.09 | 1 | 285 | .149 |
| Decision mode | .004 | [.000, .018] | 0.52 | 2 | 285 | .594 |
| Negatif affect | .005 | [.000, .027] | 1.40 | 1 | 285 | .238 |
| Positive affect | .004 | [.000, .026] | 1.25 | 1 | 285 | .264 |
| Priming \(\times\) Product | .007 | [.000, .031] | 1.89 | 1 | 285 | .170 |
| Product \(\times\) Negatif affect | .007 | [.000, .031] | 1.92 | 1 | 285 | .167 |
| Product \(\times\) Positive affect | .015 | [.000, .047] | 4.41 | 1 | 285 | .037 |
plot_model(final_price)
interactions::cat_plot(final_price, pred = Priming, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), modx.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Price Evaluation")
interactions::interact_plot(final_price, pred = positive_affect, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, modx.labels = c("A+","A++"), x.label = "Positive affect towards climate change", y.label = "Price Evaluation")
interactions::interact_plot(final_price, pred = positive_affect, modx = Product, interval = TRUE, plot.points = F, point.alpha =0.3, modx.labels = c("A+","A++"), x.label = "Positive affect towards climate change", y.label = "Price Evaluation")
\(~\)
Grand biais vers 0 avec les positive affect scale. Resultat pas très interpretable donc… Mais je laisse les deux graphs au cas ou.
\(~\) \(~\)
base_priceA = filter(base_price, Product == 1)
final_priceA = lm(Price ~ Priming*Decision_mode + Age, data = base_priceA)
apa_lm <- apa_print(anova(final_priceA))
apa_table(apa_lm$table, caption = "Anova table for Price with only A++ product.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .003 | [.000, .034] | 0.37 | 1 | 141 | .545 |
| Decision mode | .002 | [.000, .016] | 0.15 | 2 | 141 | .858 |
| Age | .004 | [.000, .039] | 0.56 | 1 | 141 | .456 |
| Priming \(\times\) Decision mode | .003 | [.000, .019] | 0.18 | 2 | 141 | .832 |
plot_model(final_priceA)
interactions::cat_plot(final_priceA, pred = Priming, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), errorbar.width = 0.4, y.label = "Price Evaluation")
base_buy = filter(base_clean, id %notin% c("13804")) #remove huge outliers, see appendix
final_buy = lm(Buy ~ Priming*Product + Age + Decision_mode + Priming*Product + Product*negatif_affect + Product*positive_affect, data = base_buy, y.label = "Probability to buy")
apa_lm <- apa_print(anova(final_buy))
apa_table(apa_lm$table, caption = "Anova table for Probability to buy.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .000 | [.000, .000] | 0.06 | 1 | 284 | .808 |
| Product | .037 | [.010, .080] | 11.03 | 1 | 284 | .001 |
| Age | .033 | [.008, .075] | 9.84 | 1 | 284 | .002 |
| Decision mode | .000 | [.000, .000] | 0.02 | 2 | 284 | .979 |
| Negatif affect | .004 | [.000, .025] | 1.14 | 1 | 284 | .287 |
| Positive affect | .002 | [.000, .019] | 0.49 | 1 | 284 | .486 |
| Priming \(\times\) Product | .000 | [.000, .001] | 0.10 | 1 | 284 | .751 |
| Product \(\times\) Negatif affect | .010 | [.000, .037] | 2.81 | 1 | 284 | .094 |
| Product \(\times\) Positive affect | .005 | [.000, .028] | 1.46 | 1 | 284 | .229 |
sjPlot::plot_model(final_buy)
interactions::cat_plot(final_buy, pred = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Probability to buy")
sjPlot::plot_model(final_buy, type = "pred", show.data = T)[3]
interactions::cat_plot(final_buy, pred = Priming, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), modx.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Probability to buy")
base_buyA = filter(base_buy, Product == 1)
final_buyA = lm(Buy ~ Priming*Decision_mode + Age , data = base_buyA)
apa_lm <- apa_print(anova(final_buyA))
apa_table(apa_lm$table, caption = "Anova table for Probability to buy only A++ product.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .001 | [.000, .021] | 0.14 | 1 | 140 | .713 |
| Decision mode | .006 | [.000, .032] | 0.40 | 2 | 140 | .672 |
| Age | .040 | [.004, .107] | 5.90 | 1 | 140 | .016 |
| Priming \(\times\) Decision mode | .006 | [.000, .034] | 0.45 | 2 | 140 | .637 |
plot_model(final_buyA)
interactions::cat_plot(final_buyA, pred = Priming, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), errorbar.width = 0.4, y.label = "Probability to buy")
\(~\) \(~\)
base_recco = base_clean #no huge outliers, see appendix
final_recco = lm(Reccomend ~ Priming*Product + Age + Decision_mode + Priming*Product + Product*negatif_affect + Product*positive_affect, data = base_recco)
apa_lm <- apa_print(anova(final_recco))
apa_table(apa_lm$table, caption = "Anova table for Reccomendation.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .002 | [.000, .020] | 0.56 | 1 | 285 | .454 |
| Product | .090 | [.044, .147] | 28.24 | 1 | 285 | < .001 |
| Age | .008 | [.000, .035] | 2.40 | 1 | 285 | .123 |
| Decision mode | .000 | [.000, .000] | 0.05 | 2 | 285 | .954 |
| Negatif affect | .014 | [.000, .045] | 4.04 | 1 | 285 | .045 |
| Positive affect | .007 | [.000, .032] | 1.98 | 1 | 285 | .160 |
| Priming \(\times\) Product | .000 | [.000, .001] | 0.00 | 1 | 285 | .944 |
| Product \(\times\) Negatif affect | .021 | [.002, .057] | 6.16 | 1 | 285 | .014 |
| Product \(\times\) Positive affect | .008 | [.000, .033] | 2.16 | 1 | 285 | .143 |
plot_model(final_recco)
interactions::cat_plot(final_recco, pred = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Reccomendation to buy")
sjPlot::plot_model(final_recco, type = "pred", show.data = T)[5]
interactions::cat_plot(final_recco, pred = Priming, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), modx.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Reccomendation to buy")
interactions::interact_plot(final_recco, pred = negatif_affect, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, modx.labels = c("A+","A++"), x.label = "Negative affect towards climate change", y.label = "Reccomendation to buy")
interactions::interact_plot(final_recco, pred = negatif_affect, modx = Product, interval = TRUE, plot.points = F, point.alpha =0.3, modx.labels = c("A+","A++"), x.label = "Negative affect towards climate change", y.label = "Reccomendation to buy")
\(~\) \(~\)
negatif_affect barely significatif and doesn’t look very informative by itself mais l’interaction mest plus interessante (avec et sans les points pour que ce soit plus digeste)
\(~\) \(~\) Example for reporting :
Product (\(F(1, 285) = 28.24\), \(p < .001\), \(\hat{\eta}^2_G = .090\), 90% CI \([.044, .147]\)) and Negative affect towards climate change (\(F(1, 285) = 4.04\), \(p = .045\), \(\hat{\eta}^2_G = .014\), 90% CI \([.000, .045]\)) affected reccomendation. However, the effect of Behav_enviro differed by Product, \(F(1, 285) = 6.16\), \(p = .014\), \(\hat{\eta}^2_G = .021\), 90% CI \([.002, .057]\).
base_reccoA = filter(base_recco, Product == 1)
final_reccoA = lm(Reccomend ~ Priming*Decision_mode + Age, data = base_reccoA)
apa_lm <- apa_print(anova(final_reccoA))
apa_table(apa_lm$table, caption = "Anova table for Reccomendation to buy only A++ product.")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .007 | [.000, .048] | 1.03 | 1 | 141 | .312 |
| Decision mode | .005 | [.000, .031] | 0.39 | 2 | 141 | .680 |
| Age | .007 | [.000, .047] | 1.01 | 1 | 141 | .316 |
| Priming \(\times\) Decision mode | .009 | [.000, .042] | 0.67 | 2 | 141 | .513 |
plot_model(final_reccoA)
interactions::cat_plot(final_reccoA, pred = Priming, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), errorbar.width = 0.4, y.label = "Reccomendation to buy")
\(~\) \(~\)
base_eco = filter(base_clean, id %notin% c("13594", "11200", "13804", "12529", "11017")) #remove huge outliers, see appendix
final_eco = lm(Ecolo_Eval ~ Priming*Product + Age + Decision_mode + Priming*Product + Product*negatif_affect + Product*positive_affect, data = base_eco)
apa_lm <- apa_print(anova(final_eco))
apa_table(apa_lm$table, caption = "Anova table for Ecological evaluation")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .000 | [.000, .000] | 0.00 | 1 | 280 | .961 |
| Product | .214 | [.148, .282] | 76.40 | 1 | 280 | < .001 |
| Age | .000 | [.000, .007] | 0.03 | 1 | 280 | .868 |
| Decision mode | .015 | [.000, .043] | 2.19 | 2 | 280 | .113 |
| Negatif affect | .000 | [.000, .000] | 0.00 | 1 | 280 | .954 |
| Positive affect | .001 | [.000, .013] | 0.14 | 1 | 280 | .707 |
| Priming \(\times\) Product | .003 | [.000, .022] | 0.76 | 1 | 280 | .384 |
| Product \(\times\) Negatif affect | .004 | [.000, .026] | 1.23 | 1 | 280 | .267 |
| Product \(\times\) Positive affect | .003 | [.000, .022] | 0.73 | 1 | 280 | .393 |
plot_model(final_eco)
interactions::cat_plot(final_eco, pred = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Ecological Evaluation")
interactions::cat_plot(final_eco, pred = Priming, modx = Product, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), modx.labels = c("A+","A++"), errorbar.width = 0.4, y.label = "Ecological Evaluation")
base_ecoA = filter(base_eco, Product == 1)
final_ecoA = lm(Ecolo_Eval ~ Priming*Decision_mode + Age, data = base_ecoA)
apa_lm <- apa_print(anova(final_ecoA))
apa_table(apa_lm$table, caption = "Anova table for Ecological evaluation only A++ product")
| Effect | \(\hat{\eta}^2_G\) | 90% CI | \(F\) | \(\mathit{df}\) | \(\mathit{df}_{\mathrm{res}}\) | \(p\) |
|---|---|---|---|---|---|---|
| Priming | .013 | [.000, .061] | 1.79 | 1 | 138 | .183 |
| Decision mode | .019 | [.000, .063] | 1.30 | 2 | 138 | .275 |
| Age | .000 | [.000, .011] | 0.02 | 1 | 138 | .890 |
| Priming \(\times\) Decision mode | .035 | [.000, .091] | 2.51 | 2 | 138 | .085 |
plot_model(final_ecoA)
interactions::cat_plot(final_ecoA, pred = Priming, interval = TRUE, plot.points = T, point.alpha =0.3, pred.labels = c("Calculation","Affective"), errorbar.width = 0.4, y.label = "Price Evaluation")
\(~\) \(~\)
\(~\) \(~\)
\(~\) \(~\)
plotNormalHistogram(base_clean$Price, main = "Price", sub = paste("skewness =", skewness(base_clean$Price, na.rm = TRUE)))
# print("Check price density by product")
# densityPlot(base_clean$Price, g = base_clean$Product) # by product
plotNormalHistogram(base_clean$Reccomend, main = "Reccomendation", sub = paste("skewness =", round(skewness(base_clean$Reccomend, na.rm = TRUE)),2))
plotNormalHistogram(base_clean$Buy, main = "Probability to buy", sub =
paste("skewness =", round(skewness(base_clean$Buy, na.rm = TRUE),2)))
plotNormalHistogram(base_clean$Ecolo_Eval, main = "Ecologic Evaluation", sub =
paste("skewness =", round(skewness(base_clean$Ecolo_Eval, na.rm = TRUE),2)))
\(~\) \(~\)
| Stepwise Model Path Analysis of Deviance Table |
| Initial Model: |
| Price ~ Priming * Product * negatif_affect + Priming * Product * |
| positive_affect + Priming * Product * Decision_mode + Priming * |
| Product * Age |
| Final Model: |
| Price ~ Priming + Product + negatif_affect + positive_affect + |
| Priming:Product + Priming:negatif_affect + Product:negatif_affect + |
| Priming:positive_affect + Product:positive_affect + Priming:Product:negatif_affect |
| Step | Df | Deviance | Resid. Df | Resid. Dev | AIC |
|---|---|---|---|---|---|
| 272 | 1492.897 | 526.9619 | |||
|
2 | 1.2988527 | 274 | 1494.196 | 523.2193 |
|
2 | 1.1771062 | 276 | 1495.373 | 519.4524 |
|
2 | 2.1924218 | 278 | 1497.566 | 515.8860 |
|
2 | 5.3835937 | 280 | 1502.949 | 512.9482 |
|
1 | 0.0454730 | 281 | 1502.995 | 510.9572 |
|
1 | 0.1449535 | 282 | 1503.140 | 508.9857 |
|
1 | 1.2190939 | 283 | 1504.359 | 507.2257 |
|
1 | 1.8189887 | 284 | 1506.178 | 505.5834 |
|
1 | 3.3639817 | 285 | 1509.542 | 504.2438 |
\(~\) \(~\)
\(~\) \(~\)
| Stepwise Model Path Analysis of Deviance Table |
| Initial Model: |
| Buy ~ Priming * Product * negatif_affect + Priming * Product * |
| positive_affect + Priming * Product * Decision_mode + Priming * |
| Product * Age |
| Final Model: |
| Buy ~ Priming + Product + negatif_affect + Age + Priming:negatif_affect + |
| Product:negatif_affect |
| Step | Df | Deviance | Resid. Df | Resid. Dev | AIC |
|---|---|---|---|---|---|
| 272 | 599.1365 | 256.7185 | |||
|
2 | 0.5872826 | 274 | 599.7238 | 253.0085 |
|
2 | 0.8112811 | 276 | 600.5351 | 249.4086 |
|
1 | 0.4022296 | 277 | 600.9373 | 247.6068 |
|
1 | 1.5140058 | 278 | 602.4513 | 246.3516 |
|
1 | 0.0234825 | 279 | 602.4748 | 244.3632 |
|
1 | 0.9743727 | 280 | 603.4492 | 242.8415 |
|
1 | 1.9725500 | 281 | 605.4217 | 241.8075 |
|
1 | 2.1206115 | 282 | 607.5423 | 240.8425 |
|
1 | 1.6186208 | 283 | 609.1610 | 239.6300 |
|
1 | 2.2333884 | 284 | 611.3944 | 238.7133 |
|
1 | 0.8585097 | 285 | 612.2529 | 237.1286 |
|
2 | 6.9999567 | 287 | 619.2528 | 236.4936 |
|
2 | 1.1043737 | 289 | 620.3572 | 233.0210 |
\(~\) \(~\)
\(~\) \(~\)
| Stepwise Model Path Analysis of Deviance Table |
| Initial Model: |
| Reccomend ~ Priming * Product * negatif_affect + Priming * Product * |
| positive_affect + Priming * Product * Decision_mode + Priming * |
| Product * Age |
| Final Model: |
| Reccomend ~ Priming + Product + negatif_affect + positive_affect + |
| Age + Priming:Product + Product:negatif_affect + Priming:positive_affect + |
| Product:positive_affect + Priming:Age + Product:Age + Priming:Product:positive_affect + |
| Priming:Product:Age |
| Step | Df | Deviance | Resid. Df | Resid. Dev | AIC |
|---|---|---|---|---|---|
| 272 | 734.5904 | 317.0503 | |||
|
1 | 0.9862800 | 273 | 735.5767 | 315.4474 |
|
1 | 0.0347883 | 274 | 735.6115 | 313.4614 |
|
2 | 9.6592597 | 276 | 745.2707 | 313.3229 |
|
2 | 1.4877211 | 278 | 746.7585 | 309.9132 |
|
2 | 2.5010516 | 280 | 749.2595 | 306.9029 |
|
2 | 0.3121104 | 282 | 749.5716 | 303.0262 |
\(~\) \(~\)
\(~\) \(~\)
| Stepwise Model Path Analysis of Deviance Table |
| Initial Model: |
| Ecolo_Eval ~ Priming * Product * negatif_affect + Priming * Product * |
| positive_affect + Priming * Product * Decision_mode + Priming * |
| Product * Age |
| Final Model: |
| Ecolo_Eval ~ Priming + Product + negatif_affect + positive_affect + |
| Decision_mode + Age + Priming:Product + Product:negatif_affect + |
| Priming:positive_affect + Product:positive_affect + Priming:Decision_mode + |
| Product:Decision_mode + Priming:Age + Product:Age + Priming:Product:Decision_mode + |
| Priming:Product:Age |
| Step | Df | Deviance | Resid. Df | Resid. Dev | AIC |
|---|---|---|---|---|---|
| 272 | 358.0852 | 104.36179 | |||
|
1 | 0.1295729 | 273 | 358.2148 | 102.46888 |
|
1 | 0.3551773 | 274 | 358.5699 | 100.76222 |
|
1 | 0.0511166 | 275 | 358.6210 | 98.80441 |
\(~\) \(~\)
\(~\) \(~\)
report::report_packages()
## - repro (version 0.1.0; Aaron Peikert, Andreas Brandmaier and Caspar van Lissa, 2021)
## - papaja (version 0.1.0.9997; Aust et al., 2020)
## - GGally (version 2.1.2; Barret Schloerke et al., 2021)
## - tinylabels (version 0.2.2; Barth, 2021)
## - apaTables (version 2.0.8; David Stanley, 2021)
## - ggplot2 (version 3.3.5; Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.)
## - stringr (version 1.4.0; Hadley Wickham, 2019)
## - forcats (version 0.5.1; Hadley Wickham, 2021)
## - tidyr (version 1.1.4; Hadley Wickham, 2021)
## - haven (version 2.4.3; Hadley Wickham and Evan Miller, 2021)
## - readr (version 2.1.1; Hadley Wickham, Jim Hester and Jennifer Bryan, 2021)
## - dplyr (version 1.0.7; Hadley Wickham et al., 2021)
## - kableExtra (version 1.3.4; Hao Zhu, 2021)
## - matrixStats (version 0.61.0; Henrik Bengtsson, 2021)
## - car (version 3.0.12; John Fox and Sanford Weisberg, 2019)
## - carData (version 3.0.4; John Fox, Sanford Weisberg and Brad Price, 2020)
## - tibble (version 3.1.6; Kirill Müller and Hadley Wickham, 2021)
## - purrr (version 0.3.4; Lionel Henry and Hadley Wickham, 2020)
## - interactions (version 1.1.5; Long JA, 2019)
## - sjPlot (version 2.8.10; Lüdecke D, 2021)
## - moments (version 0.14; Lukasz Komsta and Frederick Novomestky, 2015)
## - R (version 4.0.1; R Core Team, 2020)
## - rcompanion (version 2.4.1; Salvatore Mangiafico, 2021)
## - gtsummary (version 1.5.0; Sjoberg DD, Whiting K, Curry M, Lavery JA, Larmarange Reproducible summary tables with the gtsummary package. The R Journal 2021;13:570–80. https://doi.org/10.32614/RJ-2021-053.)
## - MASS (version 7.3.51.6; Venables et al., 2002)
## - tidyverse (version 1.3.1; Wickham et al., 2019)